% version 1.71
% date 2022/01/03
% author gxliu

function ResultJ = BnsSeiq(j0,yT,ra1,y0)

	% calculate residual error and return as Res
	function Res = ResidualError(j,ra1,yT,y0)

		% domain of definition
		tspan = [0 100];

		% initial values, S,E,I,Q	
		% y0 = [2;10;11;14];

		% four differential equations
		function dydt = vdp45(t,y,j,ra1)

			% ra1 = (r, a1)
			r  = ra1(1);
			a1 = ra1(2);

			% j = (a2,a3,b1,b2)
			a2 = j(1);
			a3 = j(2);
			b1 = j(3);
			b2 = j(4);

			% y = (S, E, I, Q)
			S = y(1);
			E = y(2);
			I = y(3);
			Q = y(4);

			Iwarn = 10;
			% if y(3) is equal to 0 (I is equal to 0), we suppose that the government 
			% will not take action. Thus both a2 and a3 are equal to 0.
			if I <= Iwarn
				a2 = 0;
				a3 = 0;
				j(1) = 0;
				j(2) = 0;
			end

			% dydt
			dydt = [-r*( b1*E + b2*I );
				 r*( b1*E + b2*I ) - E*( a1+a2 );
				 a1*E - a2*I;
				 a3*E + a2*I];

			% dydt = [-j(1)*(  j(5)*y(2) + j(6)*y(3)  );
			% j(1)*(  j(5)*y(2) + j(6)*y(3)  )   -   (  j(2)+j(3)  )*y(2);
			% j(2)*y(2) - j(3)*y(3);
			% j(4)*y(2) + j(3)*y(3)];

		end
		% ode45
		sol = ode45(@(t,y) vdp45(t,y,j,ra1),tspan,y0);
		x = [1:length(yT)];
		% estimates values 
		yG = deval(sol, x);
		yG(1,:)
		Res = sqrt(sum((yG(1,:)' - yT).^2)) / length(yT);
	end
	
	% Resierr = optimvar("Resierr");
	% convert function to optimization expressions
	% Resiexpr = fcn2optimexpr(@ResidualError, j);
	% create an optimization problem using the optimization expression
	% prob = optimproblem('Objective', Resiexpr);
	% show
	% show(prob)
	% solve
	% ResultJ = solve(prob, j0);


	f = @(j) ResidualError(j,ra1,yT,y0);
	[ResultJ,fvalue] = fmincon(f,j0,[],[],[],[],[0,0,0,0],[1,1,1,1]);
	% [ResultJ,fvalue] = fmincon(f,j0,[],[],[],[],[0.8,0.5,0,0],[1,1,0.1,0.15]);
	% ResultJ = 0;
end
